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Abstract 

A method for the prediction of acoustic scatter from 
complex geometries is presented. The discontinuous 
Galerkin method provides a framework for the devel- 
opment of a high-order method using unstructured 
grids. The method’s compact form contributes to 
its accuracy and efficiency, and makes the method 
well suited for distributed memory parallel comput- 
ing platforms. Mesh refinement studies are presented 
to validate the expected convergence properties of the 
method, and to establish the absolute levels of a error 
one can expect at a given level of resolution. For a 
two-dimensional shear layer instability wave and for 
three-dimensional wave propagation, the method is 
demonstrated to be insensitive to mesh smoothness. 
Simulations of scatter from a two-dimensional slat con- 
figuration and a three-dimensional blended-wing-body 
demonstrate the capability of the method to efficiently 
treat realistic geometries. 

Introduction 

Despite recent advances in computational aeroa- 
coustics (CAA), the simulation of acoustic scatter 
off realistic complex geometries remains problematic. 
The high-accuracy schemes developed for CAA gen- 
erally have spatial operators with large stencils that 
require smooth meshes and are poorly suited to grids 
around realistic aircraft configurations. Even the pro- 
cess of generating a block-structured mesh without the 
smoothness required for a high-accuracy method is a 
time-consuming process often measured in weeks. Un- 
structured grids about complex geometries are more 
easily generated, and for this reason, methods using 
unstructured grids have gained favor for aerodynamic 
analyses. Elowever, they have not been utilized for 
acoustics problems because the methods are generally 
low-order and incapable of propagating waves with- 
out unacceptable levels of dissipation and dispersion. 
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Attempts to extend finite-difference and hnite- volume 
methods for unstructured grids to high-order by in- 
creasing the stencil size have introduced storage and 
robustness problems. 

The discontinuous Galerkin method 1,2 is a compact 
finite-element projection method that provides a prac- 
tical framework for the development of a high-order 
method using unstructured grids. Eligher-order accu- 
racy is obtained by representing the solution as a high- 
degree polynomial whose time evolution is governed by 
a local Galerkin projection. This approach results in a 
compact and robust method whose accuracy is insen- 
sitive to mesh smoothness. The traditional implemen- 
tation of the discontinuous Galerkin method employs 
numerical quadrature for the evaluation of the inte- 
gral projections and is prohibitively expensive. Atkins 
and Shu 3 introduced the quadrature-free formulation 
in which the integrals are evaluated a-priori and ex- 
actly for a small collection of similarity elements. The 
approach has been demonstrated to possess the ac- 
curacy required for acoustics even in cases where the 
grid is not smooth. Other issues such as non-reflecting 
boundary conditions and the treatment of non-linear 
fluxes have also been studied in earlier work . 4,5 

A major advantage of the discontinuous Galerkin 
method is that its compact form readily permits a 
heterogeneous treatment of a problem. That is, the 
element topology, the degree of approximation, even 
the choice of governing equations, can be allowed to 
vary from element to element with no loss of rigor 
in the method. To take advantage of this flexibility, 
an object-oriented C++ computer program that im- 
plements the discontinuous Galerkin method has been 
development and validated. To date, however, most of 
the applications have involved benchmark problems for 
aeroacoustics 6 with relatively simple two-dimensional 
geometries. Recent work has been aimed at adding 
and validating additionally capability that is essential 
to the aeroacoustic analysis of large complex configu- 
rations. 

This paper describes the extension of the method to 
three-dimensions, the treatment of nonuniform mean 
flows, and the efficient use of parallel computing plat- 
forms. With these new capabilities, this tool will 
enable rapid aeroacoustic analyses of realistic aircraft 
configurations. When coupled with currently available 
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grid generators and large parallel computers, the en- 
tire process of mesh generation, problem setup, and 
calculation can be performed in the time span of a few 
hours or days instead of weeks or months. 


Governing Equations 

Numerical Method 

The discontinuous Galerkin method is readily ap- 
plied to any equation of the form 


<9Q 

dt 


V • F = S 


( 1 ) 


on a domain that has been partitioned into non- 
overlapping elements that cover the domain. The 
method is defined by approximating the solution in 
each element 12 in terms of an appropriate local set of 
basis functions {&/;| !<&<#} 


N 

Q « V = ^v k b k 

k =i 


and by performing a local integral projection of the 
governing equations onto the basis set in each element. 


b k ( + V • F(V) ) cm = I b k S(V) d n 


( 2 ) 


V 1 < k < N 


The coefficients of the approximate solution v k are the 
new unknowns, and the local integral projection gen- 
erates a set of equations governing these unknowns. 
Equation (2) is solved in the weak form and, for con- 
venience, is rewritten in terms of local element coordi- 
nates (£, 77 , C) 


ct 

r — 1 t-R 


<9V 

b k —-Vb k -J~ 1 F{V) ) Jdn 


■J2 J J _1 Ff(V,W )Jds = J b k S(V)Jdtt (3) 

e = l dfl e Q, 


where E denotes the number of edges, V denotes the 
trace of V on edge e, W denotes the trace of the ap- 
proximate solution in the neighboring element on edge 


e, 


d{x,y,z) 

~ d(t,r,xy 


and J = | J | . 


The term edge will be used to refer to any segment of 
an element boundary that is shared with a neighboring 
element or with the physical boundary of the domain. 
Edges will be referred to as interior edges or boundary 
edges when it is necessary to distinguish between the 
two. 

Because each element has a distinct local approx- 
imate solution, the solution on each interior edge is 
double valued and discontinuous. An approximate 


Riemann flux F^(V,W) resolves the discontinuity 
and provides the only mechanism by which adjacent 
elements communicate. The fact that this communi- 
cation occurs in an edge integral means the solution 
in a given element V depends only on the edge trace 
of the neighboring solution W, not on the whole of 
the neighboring solution W. Also, because the ap- 
proximate solution within each element is stored as 
a function, the edge trace of the solution is obtained 
without additional approximations. 

The compactness of the discontinuous Galerkin 
method has broad implications on the accuracy, ro- 
bustness, and efficiency of the method. Boundary 
conditions are easily implemented 4 either by supplying 
the approximate Riemann flux with the exact external 
solution W = W exac t (if known) , or by reformulating 
the approximate Riemann flux in terms of the interior 
solution V and the physical boundary conditions. In 
either case, the basic formulation of a boundary ele- 
ment is no different from that of an interior element. 
This is in stark contrast to high-order finite-difference 
and Hnite- volume methods in which the standard inte- 
rior operator must be reformulated near the boundary, 
usually with a negative impact on the accuracy and ro- 
bustness of the method. 

When the basis functions are polynomials of de- 
gree p, the order of accuracy of the method has been 
proven 1 to be at least p+|. In practice, the order of 
accuracy of the method is observed in most cases to 
be p + 1. In the present work, the basis set will be the 
complete set of local polynomials of a specified degree 

p: 

{b k } = {Cv j C k for 0 < i + j + k < p} 

The Riemann flux will be approximated by a simple 
Lax-Friedrechs flux of the form 

Ff (V, W) = [F(V) + F(W) - A (W - V)] /2 

where F(V) = J _1 JF(V) • n, A is greater in 
magnitude than the eigenvalues of the Jacobian of 
[F(V) + F(W)] /2, and n is an outward pointing edge 
normal. 

Discretization 

The discrete form of equation (3) is usually obtained 
by evaluating the integrals using quadrature formu- 
las of the required order 7 which is 2 p for the volume 
integral and 2p + 1 for the edge integral. Although 
this approach is straightforward, it limits the useful- 
ness of the discontinuous Galerkin method. Tensor 
products of one-dimensional quadrature formulas can 
be used to integrate quadrilateral and hexahedral ele- 
ments to any required degree; however, the number 
of terms in the quadrature summation exceeds the 
number of unknowns by a considerable margin (e.g. 
3.5 times greater for p = 4 in three dimensions). 
Aside from Dudiner’s 8 approach in which triangles and 
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tetrahedrons are mapped into quadrilaterals and hex- 
ahedrons, there are no general quadrature formulas for 
elements such as triangles or tetrahedrons. For these 
general elements, the only available quadratures are 
those that have been computed numerically and tabu- 
lated for a limited range of p. This has restricted most 
implementations of the discontinuous Galerkin method 
to quadrilateral, hexahedral, or relatively low order tri- 
angular elements. Furthermore, the implementations 
are computationally expensive. The authors are not 
aware of any discontinuous Galerkin implementations 
using tetrahedral elements. 

The quadrature-free approach 3 was developed to cir- 
cumvent this difficulty and to allow the discontinuous 
Galerkin method to be easily and efficiently imple- 
mented on general unstructured grids to any order of 
accuracy. To implement the quadrature-free approach, 
the fluxes and sources are also written as an expansion 
in terms of the basis functions: 

M M 

F(Q)^^f J (V)6 J , S^sA' 

j = i j = i 

When F(Q) is a linear function of Q, then M = N, 
and the expansion is trivial and exact. When the flux 
is a non-linear or a linear but non-constant coefficient 
function of Q, then the degree of the flux expansion 
must be at least p + 1 and M will be greater than 
N. The same comment applies to the source term S 
except that the source expansion may be truncated to 
degree p. 

Similar treatment of the approximate Riemann flux 
is only slightly more complex due to the fact, that the 
solutions on either side of an interior edge, V and W, 
are defined in terms different coordinate systems. The 
evaluation of the edge integral is simplified by rewrit- 
ing the edge trace of V and W in terms of a lower 
dimensional coordinate system associated with 

the edge: V = and W = J2k=i™ kh- 

bk denotes a lower-dimensional basis associated with 
the edge coordinate system. This, of coarse, is just 
a coordinate 'transformation, and the coefficients v 
are easily computed from v by a linear matrix op- 
erator [v] = T [v] . The trace of the flux can be 
computed either by taking the trace of the volume 
flux, [f] = T [*J -1 </f • nj , or by recomputing the flux 
from the trace of the solution. Generally, the later is 
prefered for a linear problem. Now the approximate 
Riemann flux can be expanded in terms of b as 
M 

Ff = XXA 

k = 1 
M 

= E - A (w - v)] b k /2 

k = l 

without regard for the type of element or the orienta- 
tion of the coordinate system of the adjacent elements. 


As illustrated in figure 1 , this approach allows elements 
of different type to be freely mixed in a calculation. 



Fig. 1 Edge coordinate system ({,’7) provides a 
common interface between dissimilar elements. 


Now that the functional form of the solution, source, 
and fluxes is explicit, the integrals a.nl .analytically 
evaluated to give 


d [v/,:] 
dt 


+ A- J ~ l J 


E 


+ E B ^ [ f &] = k]. 


( 4 ) 


where 


A = M -1 


J bj'VbkdCl , 

. n 


B = M -1 



and M = 

The matrices M, A, and B, depend only on the 
shape of the similarity element and the degree of the 
solution p. Thus, the set of matrices associated with a 
particular similarity element can be precomputed and 
applied to all elements that map to it at a considerable 
savings of both storage and computation. Finally, as 
a result of symmetry in the similarity elements, the 
matrices A, B e , and T e , are sparse in a predictable 
manner that is taken advantage of readily. 

Arbitrary, flat-sided, triangular and tetrahedral el- 
ements and a subset of other element types can be 
linearly mapped into a small set of similarity elements. 
A boundary element in which one or more sides are 
curved is a case in which the physical element cannot 
be linearly mapped to a similarity element. Illustrated 
in figure 2 for a triangle, the physical element can still 
be linearly mapped to a computational element, but 
now its curved sides are described by a polynomial. 
The integrals in equation (4) can still be evaluated an- 
alytically; however, the matrices are distinct for each 
element, and in general, the matrices are full. This 
wall treatment allows the elements to be sized as re- 
quired to resolve acoustic waves and not be restricted 
by the need to resolve wall curvature. Experience thus 
far indicates that only a few such elements are required 
in most cases, so the additional storage and computa- 
tion does not impact the overall calculation. 
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Fig. 2 Mapping for curved wall element. 


p 

Maximum At /A x 


Hu scheme 

5-st.a.ge scheme 

6 stage scheme 

1 

1.70 

1.61 

1.77 

2 

0.567 

0.543 

0.572 

3 

0.287 

0.271 

0.300 

4 

0.177 

0.167 

0.185 

5 

0.122 

0.106 

0.127 

6 

0.089 

0.073 

0.930 


Table 1 Stability constraints for 5 and 6 stage 
linear Runge-Kutta schemes, and Hu's 10 low- 
dispersion Runge-Kutta scheme. 


The solution is advanced in time using a Runge- 
Kutta scheme. Earlier work used the three-stage 
third-order TVD Runge-Kutta scheme of Shu and Os- 
her. 9 A stability analysis establishing the the time step 
constraints has been previously published. 3 The low- 
dissipation and low-dispersion Runge-Kutta scheme of 
Hu 10 is used in the present work. This scheme pro- 
vides fourth-order accuracy by an alternating sequence 
of 5- and 6-st.age Runge-Kutta schemes with modi- 
fied lower-order coefficients. A stability analysis of 
this scheme applied to a one-dimensional discontin- 
uous Galerkin method, given in table 1, indicates that 
the time step constraint falls between that of 5- and 
6-st.age linear Runge-Kutta. schemes. 


Physical Modeling 

The scatter of acoustic waves is well represented by 
the linearized Euler equations. Ri two-dimensions, the 
Q, F, and S of equation (1) are given by 


Q 


f 2 


{ P \ 



/ u 

P 

, F 

1 — 

0 

u 



0 
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V 0 
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0 
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0 
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o_ 
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Q 
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(5) 


An overline has been used to denote local temporal- 
mean quantities, and subscripted values denote differ- 
entiation. p and p are the density and pressure, and u 
and v are the x and y directed velocities, respectively. 
7 is the ratio of specific heats. The equations have 
been made dimensionless using the ambient, speed of 
sound c 0 as the reference velocity. The density and 
the pressure have been normalized by p 0 and p 0 c 3 , re- 
spectively. Time t is normalized by l r /c 0 where the 
reference length scale l r is case specific. In the ab- 
sence of a. mean flow, an acoustic wave will propagate 
a. unit, in the non-dimensional distance in each unit, of 
non-dimensional time. 

When the mean flow quantities are spatially con- 
stant., the flux expansion in terms of the basis functions 
by is a. trivial operation. However, when the mean 
flow quantities are non-uniform, their spatial variation 
must, be accounted for in the flux expansion. A general 
approach is to approximate the mean flow by poly- 
nomials and then derive the flux expansion in terms 
of polynomial products. To ensure the formal order 
properties of the discontinuous Galerkin method, the 
mean flow quantities should be represented to the same 
degree as the solution (degree p). The polynomial 
products can be truncated to degree p+ 1. However, 
this general approach requires considerable storage for 
the mean flow quantities. Fortunately, in wave scat- 
ter simulations, the assumption of linearity is based 
on the premise that, the perturbation amplitude is 
small relative to its wavelength (to neglect, non-linear 
steepening), and the wavelength is small relative to 
variations in the mean flow (to neglect, feedback from 
the acoustics to the mean flow). Thus, for the pur- 
poses of simulating the scatter of acoustic waves, it. is 
sufficient, to represent, the mean flow by a. lower-order 
polynomial such as a. piecewise linear representation. 


Results 


Parallel Computations 

The program has been ported to parallel comput- 
ing platforms using MPI calls. The initial port, was 
performed on the two-dimensional version of the code. 
The port, required only minimal changes to the code 
and was performed in only a. few weeks. 11 Aside 
from the mesh partitioning step, the port, involved cre- 
ating and inserting three functions containing about. 
120 lines of code: InitPid, BeginSendRecv, and 
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EndSendRecv. The InitPid function allocates the 
send and receive buffers and initializes domain connec- 
tivity data structures. The BeginSendRecv function 
loads the send buffer and posts both the MPI sends 
and receives, and the EndSendRecv synchronizes the 
calculation. The latter two functions are simply in- 
serted into the Runge-Kutta scheme at the appropriate 
points. The three-dimensional version of the code, 
which was created sometime later, ran in parallel with 
no changes to the above three functions. The speed 
and flexibility of the port is attributed to the modu- 
lar nature of the code and the language in which it is 
written (C++). 

In the two-dimensional case, the mesh is partitioned 
using the PARMETIS 12 software package. The three- 
dimensional results are preliminary, and the mesh is 
partitioned simply by distributing the list of elements 
as they are read in from a mesh file. 

Figure 3 gives the speed-up for the SP2 and clusters 
of workstations for a small case (800 third-order trian- 
gular elements). As the domain is divided over more 
processors, the number of elements assigned to each 
processor becomes too small, and t. lift communication 
overhead becomes apparent. Also, the workstations 
were not isolated or dedicated, so the performance 
figures are well below ideal; however, four to eight 
distributed workstations still provide a valuable per- 
formance improvement . 



Fig. 3 Speed-up on SP2 and workstations for a 
domain with 800 elements. 

Figure 4 gives i lie speed-up on the Origin 2000 
(02K) using an improved version of the program that 
has been optimized for cache usage on SGI machines. 
The parallel implementation has not been modified ex- 
cept for the partitioning and initialization procedures. 
Results are given for two domain sizes using fifth-order 
triangular elements (p = 4). The superlinear speed-up 
observed on the 02K is attributed to the cache per- 
formance of the code,: As the number of processors 
is increased, the domain is broken into smaller parts 


that such that, at some point, all of the data fits in 
cache. This is clearly evident in figure 5 which shows 
the computational rate as a function of the number of 
elements on each processor. For comparison, the figure 
also shows results for a range of problem sizes run on a 
single processor. Although there is a consistent degra- 
dation in performance between the single-processor 
and multi-processor cases, the computational rate of 
the multi-processor cases are similar indicating good 
scalability. Thus, computational rate may be a better 
indicator of scalability then the usual “speedup” mea- 
sure. The scatter in the data gives an indication as 
to the quality of the load balancing provided by the 
PARMETIS mesh partitioner. 



Fig. 4 Speed-up for large domain on Origin 2000. 



Fig. 5 Computational rate gives a clear indication 
of the cache performance. 

Preliminary performance results for the three- 
dimensional code are shown in figure 6 for two large 
domains. The performance is encouraging considering 
the simplistic method by which the mesh is parti- 
tioned. 
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Fig. 6 Speed-up of three-dimensional code on Ori- 
gin 2000. 

Nonuniform Mean Flows 

To validate the code for nonuniform mean flows, the 
growth rate of a temporal instability in a shear layer 
is calculated. The mean velocity profile is given by 
M( 1 + tanh(f/))/2 where the maximum Mach number 
is M = 0.2. The computational domain is comprised 
of triangles and extends from -4.25 to 4.25 in y and 
from -10 to 10 in x. Th& triangular mesh is gener- 
ated by dividing cells of a Cartesian mesh with equal 
spacing in x and sinh stretching in y. All of the calcu- 
lations are performed with degree p = 4 polynomials 
which provides fifth-order accuracy. The spatial scales 
of the mean flow and the instability wave are similar 
so it is also necessary to approximate the mean flow 
with degree four polynomials. The eigenfunction for a 
particular instability mode, calculated by the stability 
code developed by Macaraeg et a/., 13 is used to pro- 
vide initial profiles for all of the variables. Figure 7 
shows instantaneous contours for the u velocity after 
several periods of oscillation. During the simulation, 
the pressure is sampled at ten points in the vicinity of 
the shear layer. Figure 8 shows the evolution of the 
pressure at the center of the shear layer over four pe- 
riods. The growth rate from the stability code is used 
to determine the amplitude envelope which compares 
very well with the computed solution. The average 
growth rate of the pressure over the second two pe- 
riods is compared with the result from the stability 
code in figure 9. As the mesh is refined, the solu- 
tion quickly attains the value from the stability code. 
A second mesh refinement is performed by randomly 
perturbing the initial mesh by 20%) in both the x and 
y directions at each triangle vertex. Except for the 
coarsest meshes, the solution on the perturbed mesh 
is nearly as good as that on the original mesh. 

Deployed Slat 

Computations are being performed on the scattering 



Fig. 7 Instantaneous u velocity contours for a tem- 
poral instability. 



Fig. 8 Temporal evolution of the pressure at the 
center of the shear layer for a temporal instability. 

and resonan&g characteristics of a slat geometry using 
the discontinuous Galerkin method. An example grid 
is shown in figure 10. The actual grid used in the 
simulations has 6.2 times as many elements. Acoustic 
sources at various frequencies and with different ori- 
entations have been placed at several locations around 
the slat. The region of primary interest is the trail- 
ing edge of the slat where unsteady Reynolds-averaged 
Navier-St.okes (RANS) calculations 14 have shown sig- 
nificant vortex shedding. FigurfiB 11 and 12 show the 
scattered fields for dipole sources oriented normal and 
parallel to the slat edge. A frequency of 50 kHz is used 
that corresponds to the dominant frequency observed 
in the RANS calculations. There is significant beam- 
ing and a strong dependence on the source orientation. 

These calculations can be performed at a rate of 62 
seconds/period using 6 02K processors and 32 s^6- 
onds/period on 10 02K processors. With this level of 
performance, it is practical to rapidly investigate the 
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Fig. 9 Growth rate comparison for a temporal 
instability using uniform and perturbed meshes. 

effects of variations in the acoustic source and geome- 
try. 



x 


Fig. 10 Example mesh for a slat geometry. 

Three-Dimensional Simulations 

Three dimensional mesh refinement studies are per- 
formed to verify that the code produces the expected 
convergence rate as well as to determine the absolute 
levels of error one can expect for a given resolution. 
With methods for structured grids, the latter studies 
can be performed exclusively in one-dimension with 
reasonable expectations that the same error behavior 
will be observed in two or three dimensions. When 
using unstructured grids, however, the topology of the 
element changes dramatically, and it is not at all clear 
that the studies in one dimension will carry over to 
two or three dimensions. All of the three-dimensional 
results shown here are for a fifth-order method (p = 4), 

The mesh refinement tests are performed on a unit 
cube using the linear Euler equations with no mean 
flow, and an initial solution of the form p = cos(27r(i;-|- 
y)), and u = v = p/\/2. A degree four polynomial 



Fig. 11 Scattering of a dipole source at the trailing 
edge of the slat. The source is oriented normal to 
the streamlines at the trailing edge. 




* 


‘.i' 



Fig. 12 Scattering of a dipole source at the trailing 
edge of the slat. The source is oriented parallel to 
the streamlines at the trailing edge. 


approximation to the initial solution is obtained by fit- 
ting a polynomial at a collection of points within each 
element. A tetrahedral mesh is constructed from an 
N x xN y xN z Cartesian mesh in which each cube is sub- 
divided into five tetrahedrons (N x denotes the number 
of elements in the x direction) . With this construction, 
the distance between vertices varies by a factor of \/2. 
The initial test propagates the wave for 10 periods on 
an NxNxN mesh for N = 2,4, and 8. A second test 
propagates the wave for 100 periods on an NxNx2 
mesh for N = 2,4, and 8. The reported error is the 
L'> norm of the local error in the solution sampled at 
2500 points in the domain. Figure 13(a) gives the con- 
vergence rate at times = 1, 10, and 100. Between the 
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two finer grids (M = 4 and 8), all cases converge at 
a rate of about 5.5. Figure 13(b) gives the growth in 
the error as a function of time, which is equivalent to 
distance propagated. With as little at two elements 
per wavelength the error is still less than 3% at time 
= 100 . 




Fig. 13 Propagation error as a function of (a) mesh 
size, (b) distance propagated. 

It has been claimed that the discontinuous Galerkin 
method is insensitive to mesh smoothness. This prop- 
erty has been demonstrated for two-dimensional scalar 
advection 4 and for a two dimensional shear layer ear- 
lier. It is demonstrated here for three-dimensional 
wave propagation. A Gaussian pressure pulse centered 
(x, y, z) = (0.75, 0.5, 0.5) generates a spherical acoustic 
wave within a unit cube. The mean flow is uniform but 
with (f = 0.5, V' = W = 0. A tetrahedral mesh with 
6000 elements is generated from a IOaTOaTO Cartesian 
mesh in which each cube is divided into six tetra- 
hedrons. The solution is approximated with degree 
four polynomials. Figure 14(a) shows the perturba- 
tion pressure on three planes of the box at time = 


0.25. Next, each grid point is randomly perturbed by 
20%) of the mean mesh size and the case is rerun. Fig- 
ure 14(b) shows very good agreement between the two 
solutions on the line y = z = 0.5. The uniform case is 
shown as a solid line, and the case with the randomly 
perturbed mesh is shown with symbols. Contour plots 
of the solutions are indistinguishable. 




Fig. 14 Perturbation pressure of an acoustic wave 
initiated by a Gaussian pulse, a) solution from 
uniform mesh on three planes b) solution on line 
y = z = 0.5 with smooth and randomly perturbed 
meshes. 

Finally, a computation is presented on a generic 
blended-wing-body configuration in which the mesh 
is generated by commonly available tetrahedral mesh 
generation software. Figure 15 shows the surface mesh 
and several planes and a line on which results are 
shown. The aircraft has a length of 150 ft. and a 
simi-span of 145 ft. to be representative of a configu- 
ration currently under study by NASA and industry. 
The outer boundary of the computational domain is 
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a box defined by —50 < x < 200, —200 < y < 0, 
and —50 < z < 50. The geometry and surface 



Fig. 15 Generic blended-wing-body configuration 
with planes and lines on which data is examined. 


mesh are generated using the commercial software 
MSC/Patran, and volume meshes are generated us- 
ing VGRID 15 and AFLR3. lb These volume meshers 
were developed for aerodynamic performance calcula- 
tions in which the mesh size is expected to grow away 
from the body. Both meshers have difficulty gener- 
ating meshes with uniformly sized tetrahedrons. To 
ease the problem, the aircraft surface is meshed with 
an average edge length of five, and the bounding box 
is meshed with an average edge length of seven. The 
results shown here are for a mesh with 78048 tetra- 
hedrons generated by AFLR3. Ri general, meshes 
generated by AFLR3 are more uniform; however even 
the mesh used here has volumes ranging from 2.2 to 
243. 

An initial pressure disturbance of the form 



Fig. 16 Perturbation pressure along line “A” 
shown in figure 15: (a) initial transient (b) later 
solutions scaled to show the 1/r decay. 


_ f cos 2 (7ri?,/8) for R <4 

P = { 0 for R > 4 ’ 

where R = \/((x — 140) 2 + (t/ + 10 2 -|-(c — 5) 2 ), cre- 
ates a transient wave that is propagated to a time of 
160. The calculation required approximately 48 hours 
using eight processors of an Origin 2000. Based on the 
parallel performance shown earlier, it is reasonable to 
expect that such a calculation can be performed in only 
a few hours on a, large array of 80 to 100 processors. 

Figure 16(a) shows the initial evolution of this pulse 
along line “A” shown in figure 15. The pulse should 
to decay at a rate proportional to 1/7?,. Figure 16(b) 
shows the solution along line “A” at time = 50, 100, 
and 150 with an appropriate scaling. The final two 
sample times are in close agreement which indicates 
the wave is decaying at the correct rate. Figures 17(a- 
c) show the pressure on a horizontal cut through the 
winglet, at time = 130, 140, and 150. Figures 18(a- 
c) show a similar set of solutions on a spanwise cut 


through the wing and winglet. Though qualitative, 
these figures show the propagation and reflection of 
waves in the expected manner. 

Conclusions 

Several improvements to the high-order discontinu- 
ous Galerkin program have been validated. The capa- 
bility to treat non-uniform mean flows was verified by 
comparison with an instability wave on both smooth 
and random meshes. The accuracy of the solution is 
not sensitive to the smoothness of the mesh. The com- 
pact, form of the discontinuous Galerkin method allows 
easy and efficient ports to parallel platforms. The 
parallel implementation gives superlinear speed-up for 
large problems. This is attributed to the compact 
form of the discontinuous Galerkin method which al- 
lows useful work to be overlapped with communication 
as well as cache accelerations that occur when a large 
problem is broken into smaller components. The ex- 
tension to three dimensions has been performed and 
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verified for smooth and non-smooth meshes. Solu- 
tions converge to the exact solutions at the expected 
rate, and the accuracy is not greatly effected by mesh 
smoothness. Results are presented establishing the 
absolute levels of error one can expect for a given 
resolution and distance of propagation when using 
tetrahedral elements. It is demonstrated that the 
three-dimensional fifth-order method can propagate 
a wave for 100 wavelengths with less than 3% error 
while using only two elements per wavelength. Fi- 
nally, application of the program to the analysis of a 
two-dimensional slat and a generic blended- wing-body 
configuration demonstrates that the method and pro- 
gram can rapidly provide analysis of real geometries. 
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Fig. 18 Perturbation pressure on a spanwise cut 
through the wing and winglet at time = (a) 130, 
(b) 140, (c) 150. 


Fig. 17 Perturbation pressure on a horizontal cut 
through the winglet at time = (a) 130, (b) 140, (c) 
150. 
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